#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <locale.h>

#define MAX_FEATURES 100
#define MAX_SAMPLES 10000
#define MAX_ITER 1000
#define EPS 1e-12
#define PI 3.14159265358979323846

// 数据结构
typedef struct{
    double** data;
    int n_samples;
    int n_features;
}Dataset;

// 矩阵结构
typedef struct{
    double** elements;
    int rows;
    int cols;
}Matrix;

// PCA结果结构
typedef struct{
    Matrix* components;
    Matrix* explained_variance;
    Matrix* explained_variance_ratio;
    Matrix* mean;
    Matrix* singular_values;
    int n_components;
    double retained_variance_ratio;
}PCA_Result;

// 内存分配和释放函数
Matrix* create_matrix(int rows, int cols) {
    Matrix* mat = (Matrix*)malloc(sizeof(Matrix));
    if (!mat) return NULL;

    mat->rows = rows;
    mat->cols = cols;
    mat->elements = (double**)malloc(rows * sizeof(double*));
    if (!mat->elements) {
        free(mat);
        return NULL;
    }

    for (int i = 0; i < rows; i++) {
        mat->elements[i] = (double*)calloc(cols, sizeof(double));
        if (!mat->elements[i]) {
            for (int j = 0; j < i; j++) free(mat->elements[j]);
            free(mat->elements);
            free(mat);
            return NULL;
        }
    }
    return mat;
}

void free_matrix(Matrix* mat) {
    if (mat) {
        if (mat->elements) {
            for (int i = 0; i < mat->rows; i++) {
                if (mat->elements[i]) free(mat->elements[i]);
            }
            free(mat->elements);
        }
        free(mat);
    }
}

Dataset* create_dataset(int n_samples, int n_features) {
    Dataset* dataset = (Dataset*)malloc(sizeof(Dataset));
    if (!dataset) return NULL;

    dataset->n_samples = n_samples;
    dataset->n_features = n_features;
    dataset->data = (double**)malloc(n_samples * sizeof(double*));
    if (!dataset->data) {
        free(dataset);
        return NULL;
    }

    for (int i = 0; i < n_samples; i++) {
        dataset->data[i] = (double*)malloc(n_features * sizeof(double));
        if (!dataset->data[i]) {
            for (int j = 0; j < i; j++) free(dataset->data[j]);
            free(dataset->data);
            free(dataset);
            return NULL;
        }
    }
    return dataset;
}

void free_dataset(Dataset* dataset) {
    if (dataset) {
        if (dataset->data) {
            for (int i = 0; i < dataset->n_samples; i++) {
                if (dataset->data[i]) free(dataset->data[i]);
            }
            free(dataset->data);
        }
        free(dataset);
    }
}

// 矩阵操作函数
void matrix_copy(const Matrix* src, Matrix* dest) {
    if (!src || !dest || src->rows != dest->rows || src->cols != dest->cols) return;

    for (int i = 0; i < src->rows; i++) {
        for (int j = 0; j < src->cols; j++) {
            dest->elements[i][j] = src->elements[i][j];
        }
    }
}

void matrix_transpose(const Matrix* src, Matrix* dest) {
    if (!src || !dest || src->rows != dest->cols || src->cols != dest->rows) return;

    for (int i = 0; i < src->rows; i++) {
        for (int j = 0; j < src->cols; j++) {
            dest->elements[j][i] = src->elements[i][j];
        }
    }
}

void matrix_multiply(const Matrix* A, const Matrix* B, Matrix* C) {
    if (!A || !B || !C || A->cols != B->rows || A->rows != C->rows || B->cols != C->cols) return;

    for (int i = 0; i < A->rows; i++) {
        for (int j = 0; j < B->cols; j++) {
            C->elements[i][j] = 0;
            for (int k = 0; k < A->cols; k++) {
                C->elements[i][j] += A->elements[i][k] * B->elements[k][j];
            }
        }
    }
}

// 检查文件是否存在
int file_exists(const char* filename) {
    FILE* file = fopen(filename, "r");
    if (file) {
        fclose(file);
        return 1;
    }
    return 0;
}

// 检查是否为CSV文件
int is_csv_file(const char* filename) {
    const char* dot = strrchr(filename, '.');
    if (!dot) return 0;
    return strcmp(dot, ".csv") == 0;
}

// 数据读取函数
Dataset* read_csv(const char* filename) {
    FILE* file = fopen(filename, "r");
    if (!file) {
        printf("无法打开文件: %s\n", filename);
        return NULL;
    }

    // 第一遍：计算行数和列数
    int n_samples = 0;
    int n_features = 0;
    char line[4096];

    while (fgets(line, sizeof(line), file)) {
        if (n_samples == 0) {
            // 计算特征数量
            char* token = strtok(line, ",");
            while (token) {
                n_features++;
                token = strtok(NULL, ",");
            }
        }
        n_samples++;
    }

    if (n_samples == 0 || n_features == 0) {
        printf("文件为空或格式错误: %s\n", filename);
        fclose(file);
        return NULL;
    }

    printf("检测到数据: %d 个样本, %d 个特征\n", n_samples, n_features);

    // 检查是否超过限制
    if (n_samples > MAX_SAMPLES) {
        printf("警告: 样本数超过最大限制 %d，将只读取前 %d 个样本\n", MAX_SAMPLES, MAX_SAMPLES);
        n_samples = MAX_SAMPLES;
    }

    if (n_features > MAX_FEATURES) {
        printf("警告: 特征数超过最大限制 %d，将只读取前 %d 个特征\n", MAX_FEATURES, MAX_FEATURES);
        n_features = MAX_FEATURES;
    }

    // 创建数据集
    Dataset* dataset = create_dataset(n_samples, n_features);
    if (!dataset) {
        fclose(file);
        return NULL;
    }

    // 第二遍：读取数据
    rewind(file);
    int sample_idx = 0;

    while (fgets(line, sizeof(line), file) && sample_idx < n_samples) {
        char* token = strtok(line, ",");
        int feature_idx = 0;

        while (token && feature_idx < n_features) {
            dataset->data[sample_idx][feature_idx] = atof(token);
            token = strtok(NULL, ",");
            feature_idx++;
        }
        sample_idx++;
    }

    fclose(file);
    printf("成功读取数据文件: %s\n", filename);

    // 打印一些数据用于调试
    printf("数据样例(前3个样本的前5个特征):\n");
    for (int i = 0; i < (n_samples < 3 ? n_samples : 3); i++) {
        for (int j = 0; j < (n_features < 5 ? n_features : 5); j++) {
            printf("%.2f ", dataset->data[i][j]);
        }
        printf("...\n");
    }

    return dataset;
}

// 数据预处理
void compute_mean(const Dataset* dataset, Matrix* mean) {
    if (!dataset || !mean) return;

    for (int j = 0; j < dataset->n_features; j++) {
        mean->elements[0][j] = 0;
        for (int i = 0; i < dataset->n_samples; i++) {
            mean->elements[0][j] += dataset->data[i][j];
        }
        mean->elements[0][j] /= dataset->n_samples;
    }

    // 打印均值用于调试
    printf("特征均值样例(前5个): ");
    for (int j = 0; j < (dataset->n_features < 5 ? dataset->n_features : 5); j++) {
        printf("%.4f ", mean->elements[0][j]);
    }
    printf("...\n");
}

// 协方差矩阵计算
void compute_covariance_matrix(const Dataset* dataset, Matrix* cov) {
    if (!dataset || !cov) return;

    Matrix* mean = create_matrix(1, dataset->n_features);
    compute_mean(dataset, mean);

    // 中心化数据
    Matrix* centered = create_matrix(dataset->n_samples, dataset->n_features);
    for (int i = 0; i < dataset->n_samples; i++) {
        for (int j = 0; j < dataset->n_features; j++) {
            centered->elements[i][j] = dataset->data[i][j] - mean->elements[0][j];
        }
    }

    // 计算协方差矩阵: (X^T * X) / (n-1)
    Matrix* centered_T = create_matrix(dataset->n_features, dataset->n_samples);
    matrix_transpose(centered, centered_T);

    for (int i = 0; i < dataset->n_features; i++) {
        for (int j = 0; j < dataset->n_features; j++) {
            cov->elements[i][j] = 0;
            for (int k = 0; k < dataset->n_samples; k++) {
                cov->elements[i][j] += centered_T->elements[i][k] * centered->elements[k][j];
            }
            cov->elements[i][j] /= (dataset->n_samples - 1);
        }
    }

    // 打印协方差矩阵对角线用于调试
    printf("协方差矩阵对角线(前5个): ");
    for (int i = 0; i < (dataset->n_features < 5 ? dataset->n_features : 5); i++) {
        printf("%.4f ", cov->elements[i][i]);
    }
    printf("...\n");

    free_matrix(mean);
    free_matrix(centered);
    free_matrix(centered_T);
}

// 改进的幂迭代法
void power_iteration(const Matrix* A, double* eigenvector, double* eigenvalue, int max_iter) {
    if (!A || !eigenvector || !eigenvalue) return;

    int n = A->rows;

    // 更好的随机初始化
    srand(time(NULL));
    for (int i = 0; i < n; i++) {
        eigenvector[i] = (double)(rand() % 1000 + 1) / 1000.0; // 1-1000之间的随机数
    }

    // 归一化
    double norm = 0;
    for (int i = 0; i < n; i++) norm += eigenvector[i] * eigenvector[i];
    norm = sqrt(norm);
    if (norm < EPS) norm = EPS; // 防止除零
    for (int i = 0; i < n; i++) eigenvector[i] /= norm;

    double* prev_vector = (double*)malloc(n * sizeof(double));

    for (int iter = 0; iter < max_iter; iter++) {
        // 保存前一个向量
        memcpy(prev_vector, eigenvector, n * sizeof(double));

        // 矩阵乘法: v = A * v
        double* new_vector = (double*)calloc(n, sizeof(double));
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                new_vector[i] += A->elements[i][j] * eigenvector[j];
            }
        }

        // 计算特征值 (Rayleigh商)
        *eigenvalue = 0;
        for (int i = 0; i < n; i++) {
            *eigenvalue += eigenvector[i] * new_vector[i];
        }

        // 归一化新向量
        norm = 0;
        for (int i = 0; i < n; i++) norm += new_vector[i] * new_vector[i];
        norm = sqrt(norm);
        if (norm < EPS) {
            free(new_vector);
            break; // 特征值为0，提前退出
        }
        for (int i = 0; i < n; i++) eigenvector[i] = new_vector[i] / norm;

        free(new_vector);

        // 检查收敛
        double diff = 0;
        for (int i = 0; i < n; i++) {
            diff += fabs(eigenvector[i] - prev_vector[i]);
        }
        if (diff < EPS) break;
    }

    free(prev_vector);
}

// 改进的PCA实现
PCA_Result* pca_fit_simple(const Dataset* dataset, double variance_threshold) {
    if (!dataset || variance_threshold <= 0 || variance_threshold > 1) return NULL;

    PCA_Result* pca = (PCA_Result*)malloc(sizeof(PCA_Result));
    if (!pca) return NULL;

    printf("开始PCA计算...\n");

    // 计算均值
    pca->mean = create_matrix(1, dataset->n_features);
    compute_mean(dataset, pca->mean);

    // 计算协方差矩阵
    Matrix* cov = create_matrix(dataset->n_features, dataset->n_features);
    compute_covariance_matrix(dataset, cov);

    // 计算特征值和特征向量
    pca->explained_variance = create_matrix(dataset->n_features, 1);
    pca->explained_variance_ratio = create_matrix(dataset->n_features, 1);
    pca->components = create_matrix(dataset->n_features, dataset->n_features);

    Matrix* temp_cov = create_matrix(dataset->n_features, dataset->n_features);
    matrix_copy(cov, temp_cov);

    double total_variance = 0;
    for (int i = 0; i < dataset->n_features; i++) {
        double* eigenvector = (double*)malloc(dataset->n_features * sizeof(double));
        double eigenvalue;

        power_iteration(temp_cov, eigenvector, &eigenvalue, MAX_ITER);

        pca->explained_variance->elements[i][0] = eigenvalue;
        total_variance += eigenvalue;

        // 保存特征向量
        for (int j = 0; j < dataset->n_features; j++) {
            pca->components->elements[i][j] = eigenvector[j];
        }

        // 从协方差矩阵中减去当前成分 (紧缩技巧)
        if (eigenvalue > EPS) { // 只对显著的特征值进行紧缩
            for (int k = 0; k < dataset->n_features; k++) {
                for (int l = 0; l < dataset->n_features; l++) {
                    temp_cov->elements[k][l] -= eigenvalue * eigenvector[k] * eigenvector[l];
                }
            }
        }

        free(eigenvector);
    }

    // 打印特征值用于调试
    printf("特征值样例(前5个): ");
    for (int i = 0; i < (dataset->n_features < 5 ? dataset->n_features : 5); i++) {
        printf("%.4f ", pca->explained_variance->elements[i][0]);
    }
    printf("...\n");
    printf("总方差: %.4f\n", total_variance);

    // 计算方差贡献率并自动选择维度
    double cumulative_variance = 0;
    pca->n_components = 0;
    double selected_cumulative_variance = 0.0;

    for (int i = 0; i < dataset->n_features; i++) {
        if (total_variance > EPS) {
            pca->explained_variance_ratio->elements[i][0] =
                    pca->explained_variance->elements[i][0] / total_variance;
        } else {
            pca->explained_variance_ratio->elements[i][0] = 0;
        }

        cumulative_variance += pca->explained_variance_ratio->elements[i][0];
        if (cumulative_variance >= variance_threshold && pca->n_components == 0) {
            pca->n_components = i + 1;
            selected_cumulative_variance = cumulative_variance;
        }
    }

    // 如果没有达到阈值，使用所有成分
    if (pca->n_components == 0) {
        pca->n_components = dataset->n_features;
        selected_cumulative_variance = cumulative_variance;
    }

    pca->retained_variance_ratio = selected_cumulative_variance;
    pca->singular_values = NULL;

    free_matrix(cov);
    free_matrix(temp_cov);

    printf("PCA计算完成\n");
    return pca;
}

// 改进的数据转换
Matrix* pca_transform_simple(const PCA_Result* pca, const Dataset* dataset) {
    if (!pca || !dataset) return NULL;

    printf("进行数据转换...\n");

    // 中心化数据
    Matrix* centered = create_matrix(dataset->n_samples, dataset->n_features);
    for (int i = 0; i < dataset->n_samples; i++) {
        for (int j = 0; j < dataset->n_features; j++) {
            centered->elements[i][j] = dataset->data[i][j] - pca->mean->elements[0][j];
        }
    }

    // 选择前n_components个主成分
    Matrix* components_subset = create_matrix(pca->n_components, dataset->n_features);
    for (int i = 0; i < pca->n_components; i++) {
        for (int j = 0; j < dataset->n_features; j++) {
            components_subset->elements[i][j] = pca->components->elements[i][j];
        }
    }

    // 转置主成分矩阵以便乘法
    Matrix* components_subset_T = create_matrix(dataset->n_features, pca->n_components);
    matrix_transpose(components_subset, components_subset_T);

    // 转换数据: X_centered * components_subset_T
    Matrix* transformed = create_matrix(dataset->n_samples, pca->n_components);
    matrix_multiply(centered, components_subset_T, transformed);

    // 打印转换后的数据用于调试
    printf("转换后数据样例(前3个样本):\n");
    for (int i = 0; i < (dataset->n_samples < 3 ? dataset->n_samples : 3); i++) {
        for (int j = 0; j < (pca->n_components < 5 ? pca->n_components : 5); j++) {
            printf("%.6f ", transformed->elements[i][j]);
        }
        printf("...\n");
    }

    free_matrix(centered);
    free_matrix(components_subset);
    free_matrix(components_subset_T);

    printf("数据转换完成\n");
    return transformed;
}

// 结果分析
void analyze_pca_results_simple(const PCA_Result* pca) {
    if (!pca) return;

    printf("\n=== PCA分析结果 ===\n");
    printf("自动选择的维度: %d\n", pca->n_components);
    printf("保留的方差比例: %.4f (%.2f%%)\n",
           pca->retained_variance_ratio, pca->retained_variance_ratio * 100);

    printf("\n各主成分方差贡献率:\n");
    double cumulative = 0;
    for (int i = 0; i < pca->n_components; i++) {
        cumulative += pca->explained_variance_ratio->elements[i][0];
        printf("PC%d: %.4f (累计: %.4f)\n", i+1,
               pca->explained_variance_ratio->elements[i][0], cumulative);
    }
}

// 保存结果
void save_results(const char* filename, const Matrix* transformed_data) {
    if (!filename || !transformed_data) return;

    FILE* file = fopen(filename, "w");
    if (!file) {
        printf("无法创建文件: %s\n", filename);
        return;
    }

    // 写入标题行
    for (int j = 0; j < transformed_data->cols; j++) {
        fprintf(file, "PC%d", j+1);
        if (j < transformed_data->cols - 1) fprintf(file, ",");
    }
    fprintf(file, "\n");

    // 写入数据
    for (int i = 0; i < transformed_data->rows; i++) {
        for (int j = 0; j < transformed_data->cols; j++) {
            fprintf(file, "%.6f", transformed_data->elements[i][j]);
            if (j < transformed_data->cols - 1) fprintf(file, ",");
        }
        fprintf(file, "\n");
    }

    fclose(file);
    printf("降维结果已保存到: %s\n", filename);
}

// 释放PCA结果
void free_pca_result(PCA_Result* pca) {
    if (pca) {
        free_matrix(pca->components);
        free_matrix(pca->explained_variance);
        free_matrix(pca->explained_variance_ratio);
        free_matrix(pca->mean);
        free_matrix(pca->singular_values);
        free(pca);
    }
}

int main() {
    setlocale(LC_ALL, "en_US.UTF-8");

    printf("=== PCA降维工具 ===\n");

    // 直接在代码中设置CSV文件路径和方差阈值
    const char* csv_file_path = "Initial data.csv";
    double variance_threshold = 0.95; // 保留95%方差

    printf("CSV文件路径: %s\n", csv_file_path);
    printf("方差阈值: %.2f\n", variance_threshold);

    // 检查文件是否存在
    if (!file_exists(csv_file_path)) {
        printf("错误: 文件不存在: %s\n", csv_file_path);
        printf("请修改代码中的csv_file_path为您的CSV文件路径\n");
        return -1;
    }

    // 检查是否为CSV文件
    if (!is_csv_file(csv_file_path)) {
        printf("错误: 文件不是CSV格式: %s\n", csv_file_path);
        return -1;
    }

    // 读取数据
    Dataset* dataset = read_csv(csv_file_path);
    if (!dataset) {
        printf("错误: 无法读取数据文件\n");
        return -1;
    }

    // 执行PCA
    printf("\n执行PCA，目标保留%.1f%%的方差...\n", variance_threshold * 100);
    PCA_Result* pca = pca_fit_simple(dataset, variance_threshold);

    if (!pca) {
        printf("错误: PCA执行失败\n");
        free_dataset(dataset);
        return -1;
    }

    // 分析结果
    analyze_pca_results_simple(pca);

    // 转换数据
    Matrix* transformed = pca_transform_simple(pca, dataset);

    if (transformed) {
        printf("\n降维后数据形状: %d x %d\n", transformed->rows, transformed->cols);

        // 保存结果
        save_results("pca_result.csv", transformed);

        free_matrix(transformed);
    }

    // 清理内存
    free_pca_result(pca);
    free_dataset(dataset);

    printf("\nPCA降维完成!\n");
    return 0;
}